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Abstract. Dislocations are the main carriers of the permanent deformation of crystals. For simu- 
lations of engineering applications, continuum models where material microstructures are represented 
by continuous density distributions of dislocations are preferred. It is challenging to capture in the 
continuum model the short-range dislocation interactions, which vanish after the standard averaging 
procedure from discrete dislocation models. In this study, we consider systems of parallel straight 
dislocation walls and develop continuum descriptions for the short-range interactions of dislocations by 
using asymptotic analysis. The obtained continuum short-range interaction formulas are incorporat- 
ed in the continuum model for dislocation dynamics based on a pair of dislocation density potential 
functions that represent continuous distributions of dislocations. This derived continuum model is 
able to describe the anisotropic dislocation interaction and motion. Mathematically, these short-range 
interaction terms ensure strong stability property of the continuum model that is possessed by the 
discrete dislocation dynamics model. The derived continuum model is validated by comparisons with 
the discrete dislocation dynamical simulation results. 
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1. Introduction 

The plastic deformation of crystalline materials is primarily carried out by the 
motion of a large number of atomistic line defects, i.e. dislocations. Based on the 
accumulated knowledge about the behavior of individual dislocations [14], discrete dis- 
location dynamics (DDD) models [19, 9, 35, 2, 22, 12] have been well developed for the 
study of crystal plasticity in a wide range of mechanical problems. For engineer appli- 
cations, however, DDD models are limited to samples of small size (order of microns 
or below), because of their high computational costs. Hence continuum models, where 
material microstructures are represented by continuous density distributions of disloca- 
tions resulting from the local homogenization of the discrete dislocation networks, are 
practically preferred [26, 18, 24, 23, 10, 6, 1, 3, 11, 16, 34, 38, 42, 7, 15, 17, 43, 36, 28, 
31, 20, 21, 44, 4, 41, 25]. 

In order to incorporate the orientation-dependent dislocation densities and the 
anisotropic dislocation interaction and motion in the continuum model, we have em- 
ployed a pair of dislocation density potential functions (DDPFs) to describe the dislo- 
cation distribution [34, 38, 39, 37, 43]. In this representation, the intersections of the 
contour lines (of integer multiples of the length of the Burgers vector) of the two DDPFs 
4> and %l> are the locations of the dislocations, see Sec. 3 for the model in two-dimensions 
(where dislocations are infinite straight lines). Essentially, the DDPF ip characterizes 
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the local distribution of the active slip planes and the DDPF <f> restricted on a slip plane 
describes the local dislocation distribution within that plane. As a result, the derived 
continuum dislocation dynamics model takes the form of a PDE system of two DDPFs 
( t> and if), instead of equations of the single variable of scalar dislocation density in the 
existing two-dimensional continuum models in the literature reviewed above for geo- 
metrically necessary dislocations. While previous continuum model based on DDPFs 
focused on dislocation glide within slip planes [34, 38, 42, 43], the continuum disloca- 
tion dynamics equations derived in this paper incorporate both dislocation motions of 
glide and climb. The continuum dislocation model based on DDPFs also provides a 
mathematical framework for rigorous analysis of the properties of the interaction and 
dynamics of dislocations and further incorporation of other important dislocation mech- 
anisms at the continuum level (such as the Frank-Read sources [42] and dynamics of 
dislocation dipoles [4, 41]). 

In existing dislocation-density-based models of plasticity, although the long-range 
dislocation-dislocation interactions are well-captured by direct averaging, the short- 
range interactions have to be incorporated with special treatments. This is because the 
mutual interaction force between two dislocations can grow as strong as the order of 1/r, 
where r is the dislocation spacing, which leads to a strong dependence of dislocation 
dynamics on the local discrete arrangement of dislocations and further influences the 
plastic behavior of materials. However, when a discrete dislocation network is treated by 
a dislocation continuum, such short-range interactions are averaged to zero. Therefore, 
the development of continuum modelling of dislocations highly relies on effective ways 
to capture the short-range interactions on a coarse-grained scale. 

For two-dimensional dislocation configurations where all dislocations are infinitely 
straight and mutually parallel, Groma et al. [11] developed a continuum formulation 
for the short-range dislocation interaction based on a statistical approach, and such 
statistical method was further extended by Dickel et al. [5] to identify the role played 
by dislocation dipoles in crystal plasticity. However, it has been argued by comparisons 
with discrete dislocation dynamics simulations (Roy et al. [27]) that the short-range dis- 
location interaction formulas obtained based on statistic approaches do not necessarily 
apply to deterministic distributions of dislocations. 

A class of representative two-dimensional dislocation configurations widely studied 
in the literature are distributions (pile-ups) of dislocation walls consisting of straight 
and mutually parallel dislocations (e.g. [27, 29, 32, 13, 8, 40, 30, 28]). In this scenario, 
dislocation-dislocation interactions take place in both directions that are in and normal 
to the dislocation slip planes. To the best knowledge of the authors’, most available 
analytical results employing dislocation densities were obtained for regular dislocation 
wall structures, where dislocations within each wall are vertically aligned and uniformly 
spaced in the direction normal to the dislocation slip planes. For example, in their 
comparisons with results of discrete dislocation model, Roy et al. [27] also used semi- 
continuum analysis (in which discreteness normal to the slip planes are maintained) for 
pile-up of dislocation walls. With the matched asymptotic techniques, Voskoboinikov 
and coworkers [32] calculated the discrete positions of a simple dislocation structure 
formed by one horizontal row of straight dislocations near a dislocation lock, where the 
dislocation density becomes singular under a continuum setting. Hall [13] generalized 
the approach in Ref. [32] to determine the discrete positions of the dislocation walls 
of infinite length near the grain boundaries. For regular dislocation walls, Geers and 
coworkers [8] identified five regimes for the interaction energy by a single parameter 
depending on the driving force, the horizonal and the vertical spacing between neigh- 
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bouring dislocations, and they also studied the continuum limit of the equilibrium state 
in each regime as the number of the regular walls tends to infinity. The analysis in 
Ref. [8] also suggests that a single field variable describing the dislocation density is not 
sufficient for the discrete-to-continuum transition for the configuration with dislocation 
regular walls. Zhu and Chapman [40] examined the equilibria of periodically arranged 
dipole walls, and a natural transition between dipolar configurations was found con- 
trolled by the dipole height to width ratio. By investigating the local behavior of the 
mean-field stress exerted by a row of dislocations, Schulz et al. added to the contin- 
uum system a dislocation density gradient term depending on the mesh size in their 
finite element calculations [30]. Schmitt et al. derived continuum internal stress for- 
mula for dislocation glide by homogenization of dislocation microstructures under the 
assumption that the geometrically necessary dislocations form regular walls [28] . Their 
obtained formula is similar in its form to the short-range dislocation interaction term 
obtained by Groma et al. [11] using statistical approach. 

In this paper, we first systematically examine the perturbed regular edge disloca- 
tion wall structures and derive continuum short-range interaction formulas rigorously 
from discrete dislocation dynamics model by asymptotic analysis. The derived accurate 
short-range interaction formulation for such representative deterministic dislocation dis- 
tributions, together with the available results in the literature reviewed above, is able to 
give more complete understanding of the nature of the short-range dislocation interac- 
tions for parallel dislocations with the same Burgers vector in the continuum model. In 
particular, our continuum short-range formulation is expressed by higher order deriva- 
tives of the dislocation distribution; although it is similar to the continuum formula 
derived using other approaches [11, 28], the exact expressions are different. Moreover, 
by using two field variables (two DDPFs), our continuum formulation incorporates the 
anisotropy of the short-range dislocation interactions in directions along or normal to 
the dislocation slip planes, in addition to the anisotropic dislocation motions of glide 
and climb. Such anisotropy is not included in the available continuum short-range inter- 
action formulas [11, 28], and although it was examined in Ref. [8] by regular dislocation 
walls, no continuum formulation is available in the existing literature to account for 
such dislocation interaction anisotropy for general cases. 

We then incorporate these continuum short-range interaction contributions in our 
continuum PDE model. These terms are local in the sense that they depend on the first 
and second partial derivatives of the DDPFs instead of their integrals. The full contin- 
uum force (including both the long-range and short-range continuum forces) provides a 
good approximation to the discrete dislocation dynamics model. Mathematically, these 
new terms in the continuum model serve as stabilizing terms that maintain the same 
stability properties as the discrete dislocation dynamics model. Moreover, since these 
short-range interaction terms are in the form of second order partial derivatives of the 
DDPFs (f> and 'ip, they also serve as regularization terms to the continuum long-range 
force terms that are in the form of integrals of first partial derivatives of (p and ip. 

The rest of this paper is organized as follows. In Sec. 2, we reviewed the discrete 
dislocation dynamics model, from which the continuum formulation of short-range in- 
teractions will be derived. In Sec. 3, we present the continuum framework for dislocation 
walls based on the representation of dislocation density potential functions, where the 
force on dislocations consists only of the long-range Peach-Koehler force. In Sec. 4, we 
show that without short-range interactions, the continuum long-range Peach-Koehler 
force is inconsistent with the Peach-Koehler force in the discrete dislocation model for 
many common dislocation distributions. In Sec. 5, we derive continuum expressions for 
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the dislocation short-range interactions from the discrete dislocation dynamics model. 
We focus on the dislocation configurations identified in Sec. 4 where the continuum 
long-range force fails to provide stabilizing effect compared with the discrete model. In 
Sec. 6, we present the DDPF-based continuum dislocation dynamics model that incor- 
porates both the long-range and the short range continuum forces. In Sec. 7, we show 
the new continuum model is indeed able to stabilize the perturbed dislocation structures 
as the discrete dislocation model does. In Sec. 8, numerical simulations are performed 
to validate the continuum model. 

2. Discrete dislocation dynamics model 

In this section, we briefly reviewed the discrete dislocation dynamics model, from 
which the continuum formulation of short-range interactions will be derived. We con- 
sider a system of parallel straight edge dislocations, see Fig. 2.1. In this case, the 
dislocation dynamics can be reduced to a two-dimensional spatial problem, in which 
the dislocations are points in the plane orthogonal to the direction of the dislocation 
lines, that is, parallel to the 2 -axis. The Burgers vector b is along the rr-axis. The 
locations of dislocations are denoted by the points {( x m ,y n )} for integer m and n. 


JL 

j_K=>w) 


,y„) j (x m ,y„) j _(x m+ 1 ,y„) 

l_ | ( X .J»- 1) i 
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Fig. 2.1. A system of parallel straight edge dislocations. 


The Peach-Koehler force / on a dislocation is a configurational force associated with 
the change of free energy d W due to a displacement d l of the dislocation: dW = —fdl. 
The Peach-Koehler force per unit length on the dislocation is related to the stress field 
by [14] 


( @xx @ xy & xz \ 
<J yX (Tyy <J yz 
&ZX &zy &zz J 




(2.1) 


where b is the Burger’s vector with the magnitude b, r= (0,0,1) is the dislocation line 
direction, and er is the stress tensor. The component of the Peach-Koehler force in the x 
direction is parallel to the plane containing both the Burgers vector and the dislocation 
line direction (which is the slip plane), and is the glide force. The component of the 
Peach-Koehler force in the y direction is normal to the direction of the Burgers vector 
b and the dislocation line direction, and is the climb force. From Eq.( 2.1), we have the 
glide force fg = ba xy and the climb force / c = — 6c r xx . 

Using linear elasticity theory, an edge dislocation located at the point (0,0) generates 
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the following stress field [14] 


<7xv{x,y)=<Tyx{x,y) 
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(2.3) 


where /i is the shear modulus and v is the Poisson ratio, and other stress components 
vanish. 

Therefore, for a dislocation located at ( x mo ,y no ), the glide force on it generated 


(■£ mg ^m)((^mo ) (Vtiq Un ) ) 

"■( l -") [( lm 0 - Xm ) 2 + (!/ n 0 -!/») 2 ] 2 

superposition, the total glide force acting on the dislocation located at (x mo ,y no ) is 


by another dislocation located at (x m ,y n ) is 2 J^_ v) (Xm ° ’ ■ % 
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yb 2 
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Similarly, the total climb force acting on the dislocation located at (x mo , y „ 0 ) is 


fc (Xm 0 >yno) 


yb 2 


\ ' ( Vmo ym)[3{x m Q X m ) T (j/np 2/rt) ] 
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With applied stress, the total glide and climb forces acting on the dislocation located 
at (x mo , y no ) can be written as 

fg(x mo ,2/no) = /g i.x mo ,y no ) + bcr xy , (2-6) 

fc (Xm 0 , y no ) = / c dd (X m 0 , 2/n 0 ) ~ b(T° x , (2.7) 

where cr xx and a x are the components of the applied stress tensor. 

In discrete dislocation dynamics, the local dislocation velocity v is given by the 
following mobility law in terms of the Peach-Koehler force [19, 9, 35, 2] as v = M-f, 
where M is the mobility tensor and f is the Peach-Koehler force. Following [35], the 
mobility tensor can be written as M = m g (I— n®n) +?n c n(g)n, where m g is the mobility 
constant for dislocation glide, m c is the mobility constant for dislocation climb, I is the 
identity matrix, and n is the normal direction of the slip plane. For the edge dislocation 
array being considered, n=(0,l,0) T , and the dislocation velocity is given by 



where the continuum Peach-Koehler force is f= (/ g ,/ c ,0) T . 

Note that when the line direction of all the dislocation lines is changed to r = 
(0,0, — 1), the Peach-Koehler force components in Eqs. (2.4) and (2.5) do not change 
because both the dislocation line direction and the stress change their signs. In this 
case, the total glide and climb forces with applied stress in Eqs. (2.6) and (2.7) change 
to f g (x mo ,y no ) = f£ d (x mo ,y no ) - b<r° y and f c (x mo ,y no ) = fc A {x mo ,y no ) + b<j 0 xx . 
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3. Continuum formulation for dynamics of dislocation ensembles using 
dislocation density potential functions 

We consider the system of parallel straight edge dislocations as shown in Fig. 2.1. 
The number of the dislocations in the vertical direction or horizontal direction is large 
and can be considered as infinity. As all the existing continuum dislocation dynamics 
models reviewed in the introduction, our continuum model is able to describe smoothly 
varying dislocation structures and holds in an averaged sense for general dislocation 
structures by homogenizing the discrete dislocations within some representative volumes 
centered at each point [43]. 

To represent the resulting dislocation continuum, we employ a pair of dislocation 
density potential functions (DDPFs) [43] 4>(x,y) and such that, for this two- 

dimensional problem, the intersection of the contour lines 

(j){x,y) = ib and ip(x,y)=jb , (3.1) 

i,j = 0,±1,±2,---, are the dislocation lines, see Fig. 3.1. Given a smoothly varying 
dislocation structure, the local slip planes are represented by the contour lines of the 
DDPF t/j, while the dislocation lines within a slip plane are described locally by the 
contour lines of another DDPF restricted on that plane. 



Fig. 3.1. Representation of dislocation ensembles by the dislocation density potential functions 
(DDPFs). Given a smoothly varying dislocation structure, the contour line of one DDPF xft coincide 
with the slip planes, while the dislocation lines within a slip plane are given by the contour lines of 
another DDPF (ft restricted on that plane. The local average active slip plane spacing d s \ and the local 
dislocation spacing within a slip plane d in are given by Eqs. (3.3) and (3-4), respectively. Note that in 
general Weft is not necessarily normal to Wip. 


With this continuum representation of dislocation distributions, the local dislocation 
line direction is determined from the DDPFs by 


V0 x V'tp 

iiwxv^r 


(3.2) 


The local normal direction of the dislocation slip plane is in the direction of Vip, and 
the local average active slip plane spacing is given by 


dsi = 6/l|V^||. 


(3.3) 


Using the fact that the local dislocation line direction is in the direction of V</>x V^, it 
can be calculated that the local dislocation spacing within a slip plane is 

d in = 6/(||(V^x V4>) x V V’ 1 1 / 1 1 V V’ 1 1 2 ) ■ 


(3.4) 
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(In fact, d ln = 6/length of V0 in the slip plane.) 

For the two-dimensional problem considered in this paper, the Nye dislocation den- 
sity tensor is reduced to a scalar dislocation density function p(x,y), which is the number 
of dislocations per unit area [26, 18, 33]. (In fact, here the Nye dislocation density tensor 
a = p(x,y)h®T, where b = (6,0,0) and t=(0,0,1) or (0,0, — 1).) Here we define p(x,y ) 
to be the signed dislocation density: which is positive when the dislocations are in the 
+z direction and negative when they are in the — z direction. Since the local disloca- 
tion number density in the DDPF framework is d ^ t = ll v dxyi/ , || , g jg ne( j dislocation 
density can then be written as 

P( x ’y) = jp('V<f>x V0-k), (3.5) 

where k is the unit vector in the +z direction. 

For example, consider the case when the distribution of dislocations is uniform in 
y direction (normal to the slip plane) and nonuniform in x direction (within the slip 
plane). The DDPFs that describe this dislocation distribution is (p{x, y) = cp(x) and 
ip(x,y) = jj, where D is the uniform active slip plane spacing. The dislocation density 
in this case is p{x,y) = ■ 

For dislocation dynamics problems, the DDPFs (p and ip also depend on time t and 
their evolution implicitly describes the dynamics of dislocations at the continuum level, 
which is 


f 0 t Tv- V0 = o, 
\ * + v-V0 = O, 


(3.6) 


where v= (y g ,v c ) T is the local dislocation velocity at the continuum level and is calcu- 
lated from the continuum Peach-Koehler force f =(/ g ,/ c ) T following the mobility law 
in Eq. (2.8) in the two dimensional form. Here the continuum glide force f g and the 
continuum climb force f c are 


f g = f g c + (T-k)ba° xy , (3.7) 

f c = ft-(T-k)ba° xx , (3.8) 

where /g C and /^ c are the continuum glide and climb forces due to the stress field of 
dislocations, and the second term in each equation is the force due to the applied stress. 

The leading order continuum Peach-Koehler force due to the long-range dislocation 
interaction is given below in terms of the DDPFs <f> and ip, using the dislocation density 
p calculated Eq. (3.5): 


rdc,0 


(x,y) 


pb 2 

27t(1 — v) 



{x-x^Kx-x^-jy-y^ 2 } 
\{x — X \) 2 + {y — yi ) 2 } 2 


p{xi,yi)dx\dyi, 


(3.9) 


f^°{x,y) = 


pb 2 


r+oo <.+00 


2 - 7 t (1 — v) J_ c 


{y-ydmx-xxf + jy-yi) 2 } 
[(x — x \) 2 + {y — yi ) 2 } 2 


p(xi,yi)dxidyi. 

(3.10) 


These continuum long-range forces are obtained by straightforward averaging from the 
discrete dislocation dynamics model in Eqs. (2.4) and (2.5) [26, 18, 23]. 
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While previous continuum model based on DDPFs focused on dislocation glide with- 
in slip planes [34, 38, 43], the continuum dislocation dynamics equations in Eq. (3.6) 
incorporate both dislocation motions of glide and climb. Compared with the level set 
discrete dislocation dynamics method [35] in which only the intersection of the zero level 
sets of two level set functions is meaningful, the continuum dislocation dynamics equa- 
tions of the two DDPFs hold everywhere in the simulation domain, i.e. the intersections 
of all the level set pairs of the two DDPFs are meaningful here. As all the existing 
continuum dislocation dynamics models reviewed in the introduction, our continuum 
model is able to describe smoothly varying dislocation structures and holds in an aver- 
aged sense for general dislocation structures by homogenizing the discrete dislocations 
within some representative volumes centered at each point [43]. 

As to be discussed in Sec. 4, it is essential to include in continuum Peach-Koehler 
force the contributions due to short-range dislocation interactions, whose accurate ex- 
pressions will be derived in the next few sections. 

4. Inconsistency between the continuum long-range force and the dis- 
crete dislocation model 

We observe that the continuum Peach-Koehler forces based on the long-range dis- 
location interaction in Eqs. (3.9) and (3.10) are not always consistent with the forces 
from the discrete dislocation dynamics model, especially when the long-range disloca- 
tion interaction vanishes. For example, when the distribution of dislocations is uniform 
in the y direction, the dislocation density only depends on the spatial variable x, i.e. 
p(x,y) = p(x). Substituting this density into Eq. (3.9) and using j^ripjjsyzdy = 0, we 
have 


f dc ’°(x y) = — f 

/g [ ,V) 27t(1 — v) J_ 

pb 2 


+oo /*+oo 


— oo J — oo 
+oo 


(x — x\)[(x — x \) 2 — (y — yi) 2 ] 
[(x—xi ) 2 + (y — yi ) 2 } 2 


p(xi)dxidyi 


r+oo 


27t(1 — v) J_ c 
= 0 . 


{x — x\) p(x\)dx\ 


[{x-XxY-jy-yx) 2 ] 
[(x — Xi ) 2 + (y — yi ) 2 ] 2 


dy i 


(4.1) 


Thus the continuum glide force in Eq. (3.9) vanishes for this case. 

We then calculate the glide force for this case using the discrete dislocation dynamics 
model. Since the distribution of dislocations is uniform in the y direction, the locations 
of dislocations can be written as {(x m ,y no +jD)\m,j = 0,±1,±2,- • •}, where D is the 
uniform inter-dislocation spacing in the y direction. On the dislocation located at 
{x mo ,y no ), the glide force calculated from the discrete dislocation dynamics formula in 
Eq. (2.4) is 


/g {x mo ,y no ) — 


ph 2 


+oo 


27t(1 — v) 
pirb 


EE 

m j— — oo 


E 


(•Tmo ^m)[(^mo %m)~ (j-D) ] 

[{Xmo ^Xm ) 2 + {jD ) 2 } 2 


(!-^ 2 m ^ 0 cosh(27r2=V^) — r 


(4.2) 


This glide force in general is nonzero. This disagreement shows that in the continuum 
model, in addition to the leading order contribution from the long-range dislocation 
interaction, it is essential to incorporate short-range dislocation interactions at higher 
orders in the coarse-graining process from the discrete dislocation dynamics model. 
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In this paper, we will derive continuum formulas for these short-range dislocation 
interactions. We first identify all the cases in which the long-range dislocation interac- 
tion vanishes. The long-range forces are easily calculated in the Fourier space, in which 
the force formulas in Eq. (3.9) and (3.10) become 

/ g dC ’°(fcl,fc 2 ) =47T 2 bG 1 (fc 1 ,fc 2 )p(fc 1 ,fc 2 ) = - ^ ( fc 2^fcl )2 p(fci>fc2), (4-3) 

‘2nh 2 ik 3 

/c dc ’° (fcr M) = ±pi 2 bG 2 (h , k 2 )p(h M) = - 1 _ v{k 2^ k 2 ) 2 p{^ M, (4-4) 

where / is the Fourier coefficient of / of the component e^ klx+k2V \ i is the imaginary 
unit and ki,k 2 are the wave numbers. 

(i) The long-range glide force vanishes, i.e. /g C,0 = 0. This is equivalent to 
fg C ’°(ki,k 2 )=0 for any ki and k 2 . Following Eq. (4.3), if /g C,0 (fci,£; 2 ) = 0, at least one 
of the following three conditions holds for any fixed ki,k 2 : k\ =0 but k 2 y^O; k 2 = 0 but 
k\ ^ 0; or p(ki ,k 2 )=0 if k±,k 2 ^0. Thus all the solutions of ,/g c, ° = 0 are given by 

p(x,y) = J2J2P^ k2 y {klX+k2V) 

ki k 2 

= Y,p(k 1 ,o)e ik ' x +Y,p(o,k 2 y k2y 

kiy^O k 2 ^0 

= Pi(x) + p 2 (y), (4.5) 

where pi{x) and p 2 (y) are some functions. 

(ii) The long-range climb force vanishes, i.e. /c’° = 0. This is equivalent to 
fc C ’°(k\,k 2 )=0 for any k\ and k 2 . Following Eq. (4.4), if fc C ’°{ki,k 2 ) = 0, at least one of 
the following two conditions holds for any fixed ki,k 2 : k 2 = 0 but k\ ^ 0; or p(ki ,k 2 ) = 0 
if k 2 ^ 0. Thus all the solutions of /^ c ’° = 0 are given by 

p(x,y) = Y,Y,P(kuk 2 V {klX+k2v) 

ki k 2 

= Y, P(ki,0)e ikix 

ki^O 

= P3(x), (4.6) 


where pz{x) is some function. 

In summary, the long-range glide force in the continuum model vanishes if and only 
if the dislocation density has the form p(x,y) = pi(x) + p 2 (y), and the long-range climb 
force in the continuum model vanishes if and only if the dislocation density has the 
form p(x,y) = ps{x) (which means that the dislocation distribution is uniform in the y 
direction). However, the forces calculated from the discrete dislocation dynamics model 
are not necessarily zero, see the example in Eq. (4.2). In these cases, it is essential to keep 
the next order forces that represent the short-range dislocation interaction due to the 
discreteness of dislocation distributions, in the coarse-graining process from the discrete 
dislocation dynamics model. In the next section, we examine these cases and derive 
continuum force expressions to capture such short-range interactions of dislocations. 
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5. Continuum force formulation due to short-range dislocation interac- 
tions 

In this section, we derive continuum expressions for the dislocation short-range 
interactions from the discrete dislocation dynamics model. We focus on the dislocation 
configurations identified in Sec. 4 where the continuum long-range force fails to provide 
stabilizing effect compared with the discrete model. These dislocation distributions are 
uniform either within the slip planes (in the x direction) or in the direction normal to 
the slip planes (in the y direction), i.e. , 

p = p{x) or p{y). (5.1) 

We consider the dislocation configurations that are not far from a unform distribution 
(i.e. in the linear regime of the deviations). We neglect the force due to applied stress 
in this section. 

Using the representation of DDPFs described in Sec. 3, such a perturbed uniform 
dislocation wall can be described by 



^=^y+fp, 


(5.2) 


where B is the inter-dislocation spacing in a slip plane and D is the inter-slip plane 
spacing in the uniform dislocation wall. From the formula of p in Eq. (3.5), it is easy 
to calculate that Eq. (5.1) is equivalent in the linear regime to the cases that the per- 
turbations in a DDPF </> or ip are either functions of x or y 1 i.e. 


4> = 4>(x) or 4>(y) 
ip = ip(x) or ip(y). 


(5.3) 


These dislocation configurations can be summarized into four cases as shown in Fig. 5.1. 

In Case 1, the dislocation distribution is uniform in the direction normal to the 
slip planes, but nonuniform in a slip plane. This dislocation structure can be described 
using DDPFs <p and ip as (p=(p{x), ip = by/D. In Case 2, each row of dislocations has a 
small perturbation in the direction normal to the slip planes, and the perturbations are 
uniform in the direction normal to the slip planes. This dislocation structure is given 
by (f> = bx/B , ip = by/D + ip(x), where ip[x) is some function. In Case 3, the dislocation 
distribution is uniform in any slip plane, but nonuniform in the direction normal to 
the slip planes. This dislocation structure is given by <p = bx/B , i p = ip(y). Finally, in 
Case 4, each column of dislocations has a small perturbation, and the perturbations 
are uniform for all the columns of dislocations. This dislocation structure is given by 
<p = bx/B + <p(y), ip = by/D. 

We then derive for each of these four cases a continuum formula of the short-range 
dislocation interaction force from the discrete dislocation dynamics model reviewed in 
Sec. 2. In this discrete to continuum process, we assume that L»B,D where L is 
the length unit of the continuum model. This means that there are a large number of 
dislocations contained in a unit area of the domain of the continuum model. Note that 
in this limit process, b/B and b/ D are fixed finite (small) numbers, and B and D are 
greater than a few multiples of the Burgers vector length b such that the core regions 
of different dislocations are not overlapped. 

Note that although we use linear assumption, the obtained continuum model still 
holds for configurations significantly deviated from the uninform distributions. See the 
numerical examples in Sec. 8. 
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(c) Case 3 (d) Case 4 

Fig. 5.1. Four cases of dislocation distributions with vanishing long-range dislocation interaction. 
Case 1: c/) = (f)(x), 'ijj = by/D. Case 2: f> = bx/B, \f = by/D-\-'ip(x). Case 3: 4> = bx/B, = Case 

4: (p = bx / B <p(y) , 'ifj = by/D. See the text for the description of each case. 


5.1. Case 1 

The structure of dislocations in this case is shown schematically in Fig. 5.1(a), which 
is uniform in the direction normal to the slip planes (in the y direction), but nonuniform 
in a slip plane (in the x direction). This dislocation structure is described by 




^x + j)(x), ip = 'i p(y) = ^y, 


(5.4) 


where <fi(x) is some small perturbation such that (j){x)«b and (p'(x)>0. Using E- 

q. (3.5), the dislocation density p = p{x) = yj(-g + ^4^ ), and accordingly, the continuum 
Peach-Koehler force due to the long-range dislocation interaction vanishes as shown in 
Sec. 4. We will derive a continuum formula of the short-range dislocation interaction 
force from the discrete dislocation dynamics model. 

We first consider the glide force. In this case, the discrete dislocation dynamics 
model in Eq. (2.4) gives the following expression for the glide force on the dislocation 
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located at ( x m ,y n = nD ): 


f dd/_ „ \ _ / i&2 V' (X m - Xj)[{x m - Xj ) 2 - (kD) 2 } 

J g {x m ,y n ) 2n(l_v) ^ ^ - rr,U + IfcDUl 2 


j^m k =— oo 


[{x m - Xj) 2 + (kD) 2 ] 2 


(1 — v)D' 2 ^ C0S h27r^^-l 

j^m D 




•Km—j 


(1 — v)D 2 4^ \ cosh27T-^— 2p±i. — i cosh27 r- 


We will derive a continuum expression from Eq. (5.5) in the limit of the length unit 
of the continuum model L » B 1 D and b. The continuum expression will be based on 
the DDPF <j>(x) in Eq. (5.4) such that <j)(x m ) = mb, m = 0,±1,±2,---. We then have 


Xm %m+j — j B ^ [</ ) 0(x m )]. 


Using the assumption </><<&, we can make the following Taylor expansion at x m — 
Xm+j — J U ■ 


Xm-Xm+j _ -jB I B cosh27T^j--l-2Tr^j-sinh27r^ \x(„ \ _ ~ \1 

[27r ^l__rn^_ 1 - cosh27I i§.- 1 + b (cosh2-7r^- — l) 2 mXm+j) 9{X m )\ ' 


We can then approximation the glide force in Eq. (5.5) by 
fdd/„ . \^, v^+oo B cosh27r^j-n-27r^smh27r^ rI/ r 

Jg \X m ,y n )~ (1 _„) D 2 2^j=l b (cosh2^-l) 2 

Following Eq. (5.6), we have 


{xj n —j)-\-(j[Xj n +j') 2^(x m )]. 

(5.8) 


&(Xm— j ) T &(,Xm-\-j ) 2<(j(x m ) — j-.i^Xm Xm+j X m —j) — ^ (jb) Xcjxj ) — ,o j ■ 

B £> ts <p% 

(5.9) 

Note that since we have assumed (/>’(x) > 0, x can also be considered as a function of (f>. 
Thus Eq. (5.8) can be approximated by 

dd iryb 4 (f> xx ^cosh27T^-l-27r^sinh27r^ , 2 

/g {Xm,y n )tt jz i o 1 b —J 

(1 — v)D 2 4% ^ (cosh27r^-l) 2 

nfibB x C' [cosh27r^ - 1 - 2 tt£ sinh27r^](§) 2 
= l-j,^* ^ (cosh27r^ - l) 2 ’ 

irybD ( B \ 

= i Si n il. (5-10) 


1 -v aL \D 

where the function g\ (s) is defined as 


C( [27rjssinh(27rjs) — cosh(27rjs) + l](js) 2 s 
S [cosh(27rjs) — l] 2 


In the continuum model, it would be more convenient to have a simple formula for 
the coefficient instead of the summation in Eq. (5.11). Obtaining analytical formula for 



Xiaohua Niu, Yichao Zhu, Shuyang Dai, Yang Xiang 


13 


such a summation is difficult. In the following, we will derive an approximate formula 
for it. 

First, when B/D is small, the summation in Eq. (5.11) can be considered as an 
approximation to some integral with Ax = B / D as follows 


9 1 



jB_\2 
D > 


1 ^ [27r^ sinh27r^j — cosh27r^y + 1] ( 
9 ^ (cosh27 — l) 2 

(27rxsinh27rx — cosh27rx+ l)x 2 


j^o 


r*+oo 


— lim 

tc — >-0 


0 (cosh27rx— l) 2 

(27rxsinh27rx — cosh27rx+ l)x 2 B 
(cosh27rx— l) 2 D 


dx 


B 

D 


1 / r + °° 27rx 3 sinh27rx f + °° x 2 1 B\ 

2 \J-oo (cosh27rx— l) 2 X J oo cosh27rx— 1 27 t 2 D ) 

it r 2x2 dx-^y) 

2\J_ oc cosh27rx— 1 27 t 2 D J 

1 1 B 


67T 47T 2 D 


(5.12) 


Note that in these calculations, the approximation from the summation in the first line 
to the integral in the second line is based on the trapezoidal rule and the fact that the 
integrand decays exponentially as x— »±oo. Thus by Eqs. (5.10)-(5.12), we have the 
following continuum approximation of the glide force on the dislocation 




dc 

g 


db 2 ( x 3 \i/jy\ \ 
6{1 -u)\^y\\ 2n \(j) x \ ) 


(5.13) 


Here we have used -g « \<j> x \ and -jy = \4> y \ by Eq. (5.4). 

Note that the above approximation holds when B/D is small. When B/D is large, 
all the terms in the summation in g\ are exponentially small controlled by e - ^, and 
accordingly g\ is exponentially small. On the other hand, there is an important property 
that 5i>0 always holds. Thus when B/D is large, we use £r/(67 r) to approximate g i, 
where e is some small positive constant. That is, 


9i(s) 


1 S__ 

6ir 47 r 2 5 

£ 

67 r ’ 


if 1 -^ F s>e; 

otherwise for s > 0. 


(5.14) 


Fig. 5.2 shows good match between the results from the approximation of the function 
gi(s) and its exact formula in Eq. (5.11) for different values of s. 

Using the approximations in the two regimes obtained above, we have the following 
continuum approximation of the glide force on the dislocation for all values of B/ D: 


rdc 


gb 2 




1- 


_3_|V^ 

2?r | (\) x 




(5.15) 


where the notation [-] e+ is defined as 


[h\e+ 


h, if h>e\ 
£, if h< e. 


(5.16) 
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0 0.5 1 1.5 2 2.5 3 


Fig. 5.2. Comparison of the approximation of the function gi(s) in Eq. (5.14) (the red piecewise 
linear curve) and its exact formula in Eq. (5.11) (the blue curve, calculated numerically) for different 
values of s, where £ = 0.02. 


We would like to remark that in addition to its accuracy, the form of the continuum force 
formula in Eq. (5.15) is also essential to maintain the strict stability of the evolution 
equations, see Eq. (5.20). 

Note that when the line direction of these dislocations changes to r = (0,0, — 1), we 
may have (j) x < 0, and this case can be included by modifying the continuum glide force 
in Eq. (5.15) as 


rdc 

•fg 


-sgn (<j) x ) 


f-ib 2 


6(!-OI i’v 


1- 


27T \(j) x \_ 


C~\~ 


(5.17) 


where the function sgn(s) gives the sign of s. This continuum expression does not 
depend on the sign of ip y . 

Next we derive continuum expression of the climb force for this case. On the 
dislocation at ( x m ,y n ), the climb force from the discrete dislocation dynamics model in 
Eq. (2.5) is 


fc (XrrnVn) 


yb 2 y) (0 — kD)(3(x rn —Xj)) 2 + (0 — kD) 2 ) 

2tt( 1 — v) 2^ [{x m — Xj) 2 + (0 — kD) 2 ] 2 

oo J 


(5.18) 


Thus the continuum expression of the climb force in this case is 

/ c dc = 0. 


(5.19) 


Substituting the continuum expressions of /g C and f^ c in Eqs. (5.17) and (5.19) 
into the evolution equation of <f> in (3.6), with the mobility law in Eq. (2.8), the final 
form of the evolution equation for Case 1 is 




m g /r& 2 \<f) x \ 

6 ( 1 -*') \tpy\ 


JH 

2 ^ \(f>x 


: = 0 . 


£+ 


(5.20) 


5.2. Case 2 

The structure of dislocations in this case is shown schematically in Fig. 5.1(b). Each 
row of dislocations has a small perturbation in the direction normal to the slip planes (in 
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the y direction), and the perturbations are uniform in the y direction. This dislocation 
structure is described by 


4> = jjy + '4>(x), (5.21) 

where ip( x ) is some small perturbation with ij>(x)<<b and Bb/D. The continuum 
Peach-Koehler force due to the long-range dislocation interaction vanishes as shown in 
Sec. 4. 

In the discrete model of this case, if we denote the locations of the dislocations on 
the tp = 0 row by (xj = jB,yj) for j = 0,±1,±2,---, i.e., 


—y j +ii(x j ) = 0, 


(5.22) 


the glide force on the dislocation (x TO ,y m ) using Eq. (2.4) is 
fib 2 

27r(l — v) 


rdd/„ . ^ yh 2 (xm-x^Kxm-xtf-iym-iyj+kD)) 2 } 

J g {X m ,ym)- „ , x 2^ 


., , [( x rn -x j ) 2 + {y rn -{y j + kD )) 2 ] 2 

jy£mk=—oo J J 

ir fib 2 x ^ (x m — Xj)[cosh2n Xm 0 X:i cos27t Vm ^ Vi — 1] 


(1 -v)D’ 2 


E 

j^m 


(cosh27T ^j^ 2 — cos27t Vm D V3 ) 2 


Tflb 2 


(1 -v)D 2 

^ ^ ( Xm — Xj ) |^cosh2-7r I - ( cosii 2tt ' rr< j t ' J f 2) ('I — cos2tt Vrn ^ V;i j ] 


(5.23) 


Here we have summed up the contributions from each column first. When j = m, the 
glide force on the dislocation (x m ,y m ) imposed by the vertical dislocation array con- 
taining this dislocation itself is zero. The last approximation is obtained by Taylor 
expansions using the fact that cosh 27 r^ r; ^ 2 Li — 1 >> 1 — cos2 f or j ^ TO; which is 
due to Xj=jB and yj«D and B. The latter can be derived from the assumption 
ip{x) «b and Bb/D and the definition of yj in Eq. (5.22). 

Next we derive a continuum expression from the summation in Eq. (5.23) when 
B,D « L, the length unit of the continuum model. As in Eq. (5.5) in Case 1, the 
summation in Eq. (5.23) can be written in a symmetric way as 

f^( X y ^ 

J g (l_v)D 2 

(xm—Xm+j) ^cosh27r m _i_^cosh27 r m D m+ ^ ^1— cos27r Vm ^ j 

^cosh 27 r m — 1 ^ 

(Xm—Xm—j ) |^cosh27r— — — — 1— ^cosh27r— — — +2^ ^1— COs27T-^ — 

H (c OS h2-K^-^~x~iy 

(5.24) 

Using Xj =jB, Eq. (5.22), and the assumption ijj «D, we can calculation as in Case 
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1 that 


f**, x y , W & jB cosh2^+2 
/g 1 j (1 - *)U 2 (cosh2^ - l)^ 

• [V> (a: m+ j ) - rp{x m -j )] [ip ( x m + j ) - {x m ) + {x m -j )] 

Aiitt 3 D 2 - 7 P^(jB\ 4 cosh27r^ + 2 


z — ^xx^x' 

1 — 1/ ' 


/ (cosh27r^ — l) 2 


=c>0P) 


Note that we only keep linear terms of the small perturbation xp. Thus the continuum 
expression of the glide force in this case is 

/ g dc -0. (5.26) 

Next we will derive a continuum expression of the climb force in this case. The 
discrete expression given by Eq. (2.5) is 


,dd/ , n Vb 2 \ - (y m ~ (yj + kD))[ 3(x m -Xj) 2 + {y m - (y : +kD)) 2 } 
c m, m 2n(l-v) J-* [( x m -Xj) 2 -\-(y m -(yj + kD )) 2 ] 2 


jz^m k=— oo 


sin 2 TT^^i 


2(i - v)D ^ ;cosh2rr r, " ; , ' - cos27r y " 1 z / J ) 2 

. | C o S h27r^^ - cos2^^^ + 27r ^_^ S inh2^^^ J 

7r/i6 2 ^ (y m -yj)(cosh27r ;rm ~ Xj -l + 27r a;m ~ aj sinh27r Xm ~ Xj ) 


(l-^)D 2 


(cosh27r a ' m j:) Xj — l) 2 


Using the same method as before, Eq. (5.27) can be approximated by 
j. d d/_ .. 7T ybD 2 ~ ^ [cosh27T^-l + 27r^sinh27r^](^) 2 5 

fc { x m->ym) ~ /-. \ 0 4*xx / v i o jB i \o ' n • 

(l-i/jB ^ (cosh27r^-l) 2 D 


Further using IVvl = Tp |<(>x| = -g-, and taking into consider the dislocations in the oppo- 
site direction (i.e. ip y <0), as in Case 1, we have 


/ c dc = sgn(^) 


7T yb 2 \(t> x \ ( \ip. 


(1 — v)\lpy\ 2 Vl^ 


/. • j .c j f B\ v'' - ! -00 [cosh 2-7 t^?-— 1+271-^ sinh27r^](^f-) 2 g c u 

where function g<i is denned as g 2 {^) = l^ J= i 1 (cosh 2 - 7 r^-i) 2 — * d- 

stituting the obtained f^ c and f^ c into the evolution equation of in Eq. (3.6), with 
the mobility law in Eq. (2.8), the evolution equation of dislocations for this case can 

be written as ipt + S' 2 ( j|rj) ^ xx = ^ eas ^ to see 3 2 ( s ) > 0 for s > 0. 
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This means that the obtained evolution equation is not wcllposed. In order to obtain 
a wellposed equation, we can keep higher order derivative terms in the continuum ap- 
proximation, which will make the equation very complicated. Alternatively, we simply 
choose a simple regularization term of second order to ensure the wellposedness of the 
continuum model, which leads to the following evolution equation for Case 2: 

s £ipxx = 0, (5.30) 

6(1 — z') 

where £>0 is the same small parameter as that in Eq. (5.14). 

5.3. Case 3 

The structure of dislocations in this case is shown schematically in Fig.5.1(c), which 
is uniform in each slip plane (in the x direction), but slip planes of these dislocations 
are nonuniform (in the y direction). This dislocation structure is described by 

fp = j^y + '4 > {y), (5-31) 

where ip(y) is some small perturbation such that ip'(y) > 0. 

Using Eq. (3.5), the dislocation density in this case is 


Based on the conclusions in Sec. 4 (Eqs. (4.5) and (4.6)), the continuum long-range glide 
force vanishes, whereas the continuum long-range climb force does not. Therefore, in 
this case, the integral expression in Eq. (3.10) is able to give a nonvanishing leading order 
continuum approximation for the climb force, and we only need to derive a continuum 
formula for the short-range glide force. 

Using the discrete model in Eq. (2.4), the glide force on the dislocation located at 
(. x m =mB,y n ) in this case is 


fx d (x m ,y n 


y- y^ —jB[{jB) 2 — (y n — yk) 2 } 
M 1 - v) iU B ) 2 + (Vn - l Ik) 2 } 2 


(5.33) 


This means that the glide force in this case indeed vanishes. Therefore, in this case, 

/ g dc = 0. (5.34) 

Remark: In this case, we have shown that the integral expression in Eq. (3.10) is 
able to give a nonvanishing leading order continuum approximation for the climb force. 
It is interesting to note that this integral expression with the dislocation density p in 
Eq. (5.32) in this case can be further simplified to a local expression: /^ c, ° = (^^(y), 

if the perturbation if) goes to zero at infinity. 

5.4. Case 4 

The structure of dislocations in this case is shown schematically in Fig. 5.1(d). 
Each column of dislocations has a small perturbation in their own slip planes (in the 
x direction), and the perturbations are uniform in the x direction. This dislocation 
structure is described by 
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where <f>(y) is some small perturbation with cf>(y)«b and Db/B. The continuum 
Peach-Koehler force due to the long-range dislocation interaction vanishes as shown in 
Sec. 4. 

In the discrete model of this case, we denote the locations of the dislocations on the 
<p = 0 column by ( x k ,y k = kD ) for £; = 0,±1,±2,---, i.e., 

^x k + (f>(y k ) = 0- (5.36) 

The glide force on the dislocation (. x n ,y n ) using Eq. (2.4) is 


A t&2 V' (a '■n-{xk+jB))[(x n -{x k +jB)) 2 -{y n -y k ) 2 ] 

Jr \ x n,yn) — ■. / , / y 


k^nj =— oo 


[(x n ~(x k +jB)) 2 + (y n -y k ) 2 ] 2 


27t(1 — v) 
yb 2 

2(1 — v)B (cosh27 — cos27 r 1 "^ ) 2 


-E- 


sm2TT Xrl B Xk 




cosh27r^-^-cos27r^-^-27r^-^sinh27r yTt Vk 


yb 2 


B B B B 

s\\\2'K Xn -j 3 Xk (cosh27r yn g i ' fc — 1 — 2ir Vn B Vk shih2ir Vn B Vk ) 


2(1 -v)B 


E 

k^n 


(cosh2^ 


Vn-Vk 


-l ) 2 


(5.37) 


Here we have summed up the contributions from each row first. When k = n, the glide 
force on the dislocation (x n ,y n ) imposed by the row of dislocations containing this 
dislocation itself is zero. The last approximation is obtained by Taylor expansions using 
the fact that cosh27r i/n g 2/ *’ — 1 >> 1 — cos27r a: ")) a:fc for k^n, which is due to y k = kD 
and x k «B and D. The latter can be derived from the assumption (j>(y)«b and 
Db/B and the definition of x k in Eq. (5.36). The relative error of this approximation is 
0(max|aife|/ll) 2 . 

Following Eq. (5.36), we have the Taylor expansion that 


•&k %n — 


B ~ B ~ 

j4>(yk) + -^<t>(.yn) 

B ~ B ~ 

j4>v(yn)(yk - y n ) - ^ 7 4> vv (yn)(yk - y n f + 0((y k 


yn) 3 ). 


Using Eqs. (5.37) and (5.38) and y k = kD, we have 


(5.38) 


cosh27 :*§- - 1 - 2tt^ sinh27r^ 

(cosh27r^ — l) 2 

7i -yb 2n( (cosh27r^-l-27r^sinh27r^)(fcL)) 2 

“(1 ~v)B^ {yn) f^ (cosh27r*£-l) 2 


Using the definition of the function gi in Eq. (5.11) and the approximation in 
Eq. (5.14), we have the continuum approximation 


/« 


dc 


7 rybB 2 f D\ ~ 

(1 -v)D 91 \B ) ^ vv 


yb 2 \l/)y\ 

6(l-v)\<j) x \ 2 


1 - 


3 | <j>x 


2tt \ip y \\ 


£+ 


(5.40) 
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Here we have used <f) yy = (f> yy . 

As in Case 1, when the line direction of these dislocations changes to r = (0,0, — 1), 
we may have 4> x < 0, and this case can be included by modifying the continuum glide 
force in Eq. (5.40) as 


rdc 

•'g 


-sgn (<f> x ) 


^b 2 \i/J v \ 

6(1-^)! <t>x\ 2 



3 \4> x 


27 r | ip y 


l/ yy 




(5.41) 


This continuum expression does not depend on the sign of ifiy. 

As in the previous cases, we also calculate the continuum approximation of the 
climb force in this case from the discrete model in Eq. (2.5), and the result is 


,dd/, , \ yb 2 ST' V' (.Vn - Vk)[ ~ ( x k + j B)) 2 + (y n - y k ) 2 ] 

c Xn,Vn 27t(1 — v) [( x n -(x k +jB)) 2 + (y n -y k ) 2 ] 2 

fib 2 1 

2(1 — v)B (cosh27r y "g yfc —cos2n Xn g Xk ) 2 

•[— 2n Vr, g yk (cosh27r yn ~ yfc cos27r :r,> ~ a:fc — 1) 

+2sinh2ir Vr, g Vk (cosh27 j Vrl ~^ k — cos27r )] 

= 0(4>'(y n )<t>"(yn )) 

«0. (5.42) 

Again, we have used the fact that cosh27r y "~ yfc — 1 >> 1 — cos27r x " j? j;>f for k^n, to 
obtain the expansions. Therefore, in this case, 

/ c dc = 0. (5.43) 


Substituting Eqs. (5.41) and (5.43) into Eq. (3.6), we have the following evolution 
equation for this case: 




m g fib 2 | il) y \ L 

6(1 — 77 ) \<j) x \ [ 


3 \U 

2ir \i/j y \_ 



(5.44) 


6. Continuum dislocation dynamics model incorporating short-range in- 
teractions 

In this section, we present the continuum dislocation dynamics model that incor- 
porates the short range dislocation interactions discussed in the previous section. 

6.1. The continuum dislocation dynamics model based on DDPFs 

We have shown in Sec. 4 that a continuum model with only the long-range Peach- 
Koehler force is not always able to capture the behaviors of discrete dislocation dy- 
namics. It will be shown in Sec. 7 that such inconsistency leads to insufficiency in the 
stabilizing effect of the continuum model compared with the discrete dislocation dynam- 
ics model. As a result, in numerical simulations using such a continuum model, there 
is no effective mechanism to eliminate some numerical oscillations generated during the 
simulations. 

In Sec. 3, we have presented the framework of our DDPF-based continuum dis- 
location dynamics model, see Eq. (3.6). We incorporate into our continuum model 
the continuum short-range forces obtained in the previous section for the cases where 
the continuum long-range Peach-Koehler force vanishes. With these short-range terms 
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and including the contributions of the applied stress field, the continuum dislocation 
dynamics equations in Eq. (3.6) become 


>+v-V^=^^ : 


XX 1 


nigiib 2 \ipy\ L 

6(.l— i/) |0x| 


3 <t>x 
2n tp y 


£+ 


4>yy, 


( 6 . 1 ) 


where 


v= (tig, ti c ), 

tig = ^-(r-k) G\ * (y<j ) x V^-fe) +TO g (r-k)6cr! 


( 6 . 2 ) 


xyi 


= 2n(l^v) ( T ' k ) f-tt f-oo ^ WtfaM) X Vt/l(ail,t/i)-fc] d Xl d yi 


27r(l— u ) v’ sj — oo 

+m g {r -k)ba° xy , 


2 +(y-ni) 2 ] 


(6.3) 
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( 6 . 8 ) 


Under the assumptions that the length of the Burgers vector b«L , where L is the 
unit length of the continuum model, and the average dislocation spacing B ~ D «L 1 
it is easy to find that the ratio of the second order partial derivative terms vs the long- 
range terms v - Vtfi and v- Vi/’ in Eq. (6.1) is 0(b/L) « 1. Here we have used the fact 
that (j) x = 0{b/B), (f> xx = O (b / (B L)) , and similar orders for other partial derivatives of 
( t> and ip. 

Recall that continuum short-range interaction terms provide good approximations 
to the discrete dislocation model when the continuum long-range force vanishes for some 
non-trivial perturbed dislocation walls. For a general dislocation distribution described 
by the continuum model, the full continuum force (including both the long-range and 
short-range continuum forces) still provides a good approximation to the discrete dis- 
location dynamics model under the assumption that a point in the continuum model 
corresponds to one of these dislocation microstructures of perturbed regular dislocation 
walls, which is a common technique for the coarse-graining from micro- or meso-scopic 
models to continuum models. Mathematically, these short-range terms in the continu- 
um model serve as stabilizing terms that maintain the same stability properties as the 
discrete dislocation dynamics model, as will be shown in Sec. 7 

Recall also that the main advantage of continuum model based on DDPFs [34, 38, 
42, 43] is being able to describe the orientation-dependent dislocation densities of curved 
dislocations. The dislocation glide within its slip plane due to the long-range Peach- 
Koehler force is regularized by the local curvature term due to line tension effect. In 
the continuum dynamics equations in Eq. (6.1) obtained in this paper, the short-range 
interaction terms are in the form of second partial derivatives of the DDPFs and are 
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able to provide regularization in the cross-section of the dislocations for both glide and 
climb. Combining these two regularization effects, we expect to have a full well-posed 
continuum dislocation dynamics model based on DDPFs. Moreover, the use of two 
DDPFs 0 and 0 in the continuum dislocation dynamics model enables the study of the 
anisotropic behaviors of dislocation ensembles within and out of their slip planes. These 
will be further explored in the future work. 

6.2. Continuum model for dislocation glide In this subsection, we consider 
the dynamics of dislocations only by their glide. In this case, we assume the average 
inter-slip plane distance is D [43], that is, %j)(x,y) = -^y is always fixed. Applying our 
continuum model in Eq. (6.1) to this case, i.e. , the dislocations only move in the x 
direction. In this case, Eq. (6.1) becomes 


0t + ^ | (\> x | Gi * (j) x + m g ba° y \ 0 : 
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In this case, the continuum model in Eq. (6.9) can be written as: 


4>t+v g <j) x = 0, (6-11) 

where the total glide velocity v g = 7n g f e , the continuum total glide force / g = /g C + 
sgn {4>x)bcr xy as given by Eq. (3.7), and the continuum force due to interactions between 
dislocations 
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including both the long-range interaction force (the first term) and the short-range 
interaction forces (the remaining two terms) on the dislocations. 

When the dislocation distribution is uniform in the y direction, which is Case 1 in 
Sec. 5, Eq. (6.9) reduces to 


<j) t +m g ba° xy \(l) x \ 


m g yb 2 D 
6(1 — v) 


\<t>x\ 



3 b 

2ttD\(I> x \ 



(6.13) 


In this case, the continuum total force in Eq. (6.12) reduces to Eq. (5.17). 

6.3. Comparison with scalar dislocation density based continuum models 


In this subsection, we examine the evolution of the signed dislocation density p 
defined Eq. (3.5) in terms of the DDPFs 0 and 0. 

We first consider the continuum model of 0 and 0 in the form of Eq. (3.6). From 
Eqs. (3.5) and (3.6), we can calculate that 


Pt + V • (pv) = 0, 


(6.14) 
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where v= ( v g ,v c ) T is the dislocation velocity. In fact, 

Pt 'y^ipfx^fy '4 J x4 > y')t 

= Ty (f’xt'lfy T t^xt^y ^x^yt) 

= -^fyxlpy + {-V ■ Vlp)y4> X - (-V- Vlf)) x 0y - ( -V -V(/))y1p X } 

= T2{(- V 10x'lfiy + Vl'lt’x</>y)x + (-V2<fix‘ipy + V2l/>x<t'y)y} 

= -V-(pv). (6.15) 

In most of the continuum dislocation dynamics models in the literature, the evolution 
equation is written in the conservative form in Eq. (6.14). Here we only consider the 
geometrically necessary dislocations. When only the long-range Peacli-Koehler force is 
considered, the dislocation velocity v is expressed by the mobility law in Eq. (2.8) and 
the long-range force f=(/ g ,/ c ) T in Eqs. (3.9) and (3.10) in terms of p. These form a 
closed evolution equation for the dislocation density p. 

However, the modified continuum dislocation dynamics models incorporated with 
short-range interaction terms in Eq. (6.1) in general is not able to be described fully 
by the evolution of p. The reason is that in our continuum model incorporates the 
anisotropy of dislocation structure and motion within and out of the slip planes, where- 
as the single scalar dislocation density p is only able to describe isotropy dislocation 
structure and motion. When we only consider the glide motion of dislocations as in 
Sec. 6.2, following Eq. (3.5), the dislocation density is p= (V^ x Vi/’- k)/b 2 = and 

Eq. (6.11) can be written as 


Pt + {pv g)z = 0, 


(6.16) 


where v g = m g f g , /g = f g + sgn(p)6<r° , and 
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Although Eq. (6.16) is in a conservative form of the dislocation density p, the continuum 
total glide force in Eq. (6.17) also depends on <f> yy , which cannot be simply expressed 
in terms of p. Especially, for the dislocation structure of Case 4 shown in Fig. 5.1(d), 
the dislocation density p=l/(BD) : thus the representation by p alone is not able to tell 
the difference between this dislocation structure and a uniform distribution. 

For the dislocation structure of Case 1 shown in Fig. 5.1(a) (without the applied 
stress), our continuum model in Eq. (5.20) can indeed be rewritten as an evolution 
equations of the dislocation density p following p=^(f) x , which is 
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In this case, only the local short-range force is nonvanishing, which is 
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(6.18) 


(6.19) 
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In the available continuum formulas in the literature for this case using different meth- 
ods [11, 28], their local forces are proportional to p x /\p\ when only the geometrically 
necessary dislocations are consider. The corresponding term in our continuum model 
for this case in Eqs. (6.18) and (6.19) is p x D 2 , which means that for this special case, 
the isotropic dislocation density p in the denominator in the models in the literature 
should be replaced by a more accurate expression 1/D 2 where D is the average inter- 
dislocation distance normal to the slip plane. Again we can see that our model using 
two DDPFs cp and ip are able to anisotropy of dislocation structure and motion within 
and out of the slip planes, which is not able to be achieved by using the traditional 
scalar dislocation density p. 

7. Stability using the new continuum model 

In this section, we examine the stability of the uniform dislocation distributions 
using the derived continuum model in Eq. (6.1). Consider a uniform distribution of 
dislocations represented by 4>o — -gX, ipo = jjV- This uniform distribution subject to a 
small perturbation can be written as 


f <P=^x + (p(x,y,t), 


where (p(x,y,t) and ip(x,y,t) are small perturbation functions. Using Eq. (3.5), the 
dislocation density for this distribution up to linear order of the small perturbations is 
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Substituting the above (p and ip into the continuum model in Eq. (6.1) with E- 
qs. (6.2)-(6.8), the linearized evolution equations of <p(x,y,t),ip(x,y,t), written in the 
Fourier space, is 
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(7.3) 

(7.4) 


where k± and /c 2 are frequencies in the x and y directions, respectively. Here we have 
used Gi(k 1 ,k 2 ) = -i 2 ^-,) (#+%? and ^(fci.fa) = for G x {x,y) 

and G 2 {x,y) in Eqs. (6.7) and (6.8). 

The evolution of <p and ip described by Eqs. (7.3) and (7.4) is determined by the 
two eigenvalues of the coefficient matrix solved from the characteristic polynomial 
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(7.5) 
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The Fourier coefficients of the small perturbations <p and ip decay when the two eigen- 
values Ai , A 2 < 0 . 

Due to AS = RC, the characteristic polynomial in Eq. (7.5) becomes 

A 2 + (Al + ai +02 + 5 + si)A+ («j +a 2 )(S'-l-si) + Asi = 0. (7 . 6 ) 


Thus the two eigenvalues are 

x — (A+ai+a2+<S , +si)i-^/ (A+ai+a2+<S , +si) 2 — 4[(ai+a2)(*S , +si)+Asi] 

A l,2 = 2 

- (A+ai +a2+'S , +^i (A-\-ai-\-a,2 — S— si) 2 + 424 .S 

= 2 ’ 


(7.7) 
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Note that A,S,a\,a 2 ,s\ >0. By Eq. (7.8), we know that both Ai and A 2 are real. It is 
easy to conclude from Eq. (7.7) that Ai i 2<0 when Aq^O or fc 2 ^ 0 (because the term 
4 [(ai+a 2 )(S + si)+.Asi] >0 in this case), and Ai = A 2 = 0 when Aq = Aq = 0. Therefore, 

when (Aq,Aq) 7 ^ (0,0), <p and ip always decay and the uniform distribution of dislocations 
is stable using the derived continuum model in Eq. (6.1). 

This stability result provides a basis for wellposedness of the continuum model in 
Eq. (6.1) as well as stability of numerical solutions for it. These topics will be further 
explored in the future work. When only the continuum long-range Peach-Koehlcr force 
is considered, i.e. , the second partial derivative terms on the right-hand side of the PDE 
system in Eq. (6.1) vanish, the linearized equations for the small perturbations <p and 
ip in the Fourier space are 
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(7.9) 

(7.10) 


Same as the discussion in Sec. 4, when Aq = 0 or Aq=0, there is no stabilizing force 
(which is the glide force) for <p ; and when Aq = 0, there is no stabilizing force (which is 
the climb force) for ip. In these cases, numerical oscillations in simulations cannot be 
stabilized without the second order partial derivative terms. 

Recall that the second order partial derivative terms in Eq. (6.1) are based on 
the short-range interactions of dislocations. Those terms coming from the short-range 
glide forces (in the ^-equation) agree with the glide forces using the discrete dislocation 
model for uniform dislocation distributions subject to small perturbations in the glide 
direction. For the climb force, a regularization term (in the ^-equation) is added in 
addition to the stabilizing effect provided by the long-range climb force. 
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8. Numerical simulations 

In this section, we perforin numerical simulations to validate the derived continuum 
model. In addition to the nondimensionalization before simulations, we set Poisson 
ratio jz= 1/3. 

8.1. Comparisons of the continuum force with the discrete model 

We first examine the total glide force in the continuum model including both the 
long-range and short-range contributions given by Eq. (6.12) by comparisons with the 
discrete dislocation dynamics model. Recall that the continuum short-range glide force 
terms are derived from the discrete dislocation model for uniform dislocation distribu- 
tions subject to small perturbations in the glide direction. 

We first consider the dislocation distributions of Case 1 in Sec. 5, where the dis- 
location distributions are uniform in the y direction. This problem is reduced to a 
one-dinrensional problem depends only on x. 

Example 1 

Assume the dislocation distribution is described by 

4>{x)= \ ^erf(f) if-^<a;<4£ (8.1) 

W if x = 

where erf(a;) = -^= Jq e ~ u du, and V’(y) = i>2/- Periodic boundary condition is assumed 
in the x direction. We set D = 506, 13 = 305, and IV = 40. The dislocation walls are 
concentrated within the region in the center with width w. We perform simulations for 
the cases of w= 1013, w = f>B, w = B. The profiles of the DDPF <f> and the locations of 
the dislocation walls are shown in Fig. 8.1 (a), (c), and (e), and the corresponding glide 
forces calculated by the continuum model in Eq. (6.12) (which reduces to Eq. (5.17) in 
this case) and by the discrete dislocation model are plotted in Fig. 8.1 (b), (d), and (f), 
respectively. It can be seen that the results of the continuum model agree excellently 
with those of discrete model for smoothly varying (the case of w = 10B in Fig. 8.1 
(a),(b)) and even concentrated (the case of w = 5B in Fig. 8.1 (c),(d)) distributions of 
dislocation walls. For extremely concentrated distribution of dislocation walls as shown 
in Fig. 8.1 (e) with w = B, the overall continuum approximation is still very good except 
at the two ends of the concentrated distribution where the dislocation density changes 
dramatically. 

Example 2 

In this example, we examine the continuum glide force in Eq. (5.17) for different 
values of the ratio B/D for distributions of dislocation walls with uniform active slip 
plane spacing. Recall that B is the average inter-dislocation distance within a slip plane 
and D is the average slip plane spacing. For these dislocation distributions, we choose 
the DDPFs ip{y) = y and />( x) determined by the following equation 

! X = ^ + * Si "(ll)- (8.2) 

This a uniform dislocation wall distribution with perturbation in the x direction, and 
the DDPF </> can be written as 4>{x ) = ^x + (j>{x) 1 where 4>{x) is a small perturbation, see 
Fig. 8.2(a). The period of this distribution is IV = 40 dislocation walls. We fix D = 50b 
and vary the value of B. 
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Fig. 8.1. Example 1: Continuum glide force compared with the discrete model for distributions of 
dislocation walls for different values of concentration width w (defined in Eq. (8.1 ) ). The concentration 
width w = 10B in (a) and (b), w = 5B in (c) and (d), andw = B in (e) and (f). Images (a), (c), and 
(e) show the profile of (f)(x) (red curve) and locations of the dislocation walls for each value of the 
width w. The black dots on the horizontal line indicate the locations of the dislocation walls , and the 
blue dots show the corresponding values of (f in the continuum model. Images (b), (d), and (f) show 
values of the glide force / g on the dislocation walls calculated by using the continuum model (in red) 
and by using the discrete dislocation model (in blue). 


The values of the glide force calculated by the continuum model and comparisons 
with the results of the discrete dislocation model are shown in Fig. 8.2(b)-(f) for the cas- 
es of r? = 156, 406, 506, 1006, 2006. When the inter-dislocation wall distance B is smaller 
than the slip plane spacing D , as shown in Fig. 8.2(b), the continuum glide force agrees 
excellently with the force in discrete model. In this case, the glide force is significant: 
around 10~ 3 /i&, in agreement with the strong interaction between neighboring dislo- 
cation walls. When the inter-dislocation wall distance B is comparable with the slip 
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(d) (e) (f) 


Fig. 8.2. Example 2: Continuum glide force compared with that of the discrete model for distri- 
butions of dislocation walls for different values of the ratio B/D for distributions of dislocation walls 
with uniform active slip plane spacing (given by Eq. (8.2 ) ). (a) The profile of the DDPF (f> (red curve) 
and locations of the dislocation walls. The black dots on the horizontal line indicate the locations of the 
dislocation walls, and the blue dots show the corresponding values of cf in the continuum model. Images 
(b)-(f) show the continuum glide force (red dots) compared with the force calculated from the discrete 
dislocation model (blue dots) for the cases of (b) B = 15b, (c) B = A0b, (d) B = 50b, (e) B = 100b, and 
(f) B = 200b, respectively. 


plane spacing D , as shown in Fig. 8.2(c) and (d), the continuum glide force agrees well 
with the force in discrete model with small errors. In this case, the glide force becomes 
smaller: around 10 _4 /i6, which is again consistent with the weak interaction between 
neighboring dislocation walls in this case. When the inter-dislocation wall distance B 
is much greater than the slip plane spacing D, the interaction between neighboring dis- 
location walls should be negligible, which is reflected by the small values of the forces 
calculated by the continuum and the discrete models shown in Fig. 8.2(e) and (f) (at 
the order of < 10~ 6 /i6 and < 10~ lo /i6). In this sense, the continuum model still provides 
a good approximation to the discrete model in this case, although the values calculated 
by the two models are not necessarily exactly the same. The latter differences at the 
negligible orders are due to the simplification of our continuum model in Eq. (5.17) from 
its exact form in Eqs. (5.10) and (5.11) using the simplification in Eq. (5.14). 

Example 3 

In this example, we examine the continuum glide force in Eq. (6.12) for a general 
dislocation distribution. The dislocation distribution is given by 

f ^(a;,2/) = -|a ; + 0.02sin(^10x)sin(|i2j/), 

{ x,y ) = ^y + 0.02sin(|^2a;) sin (f^5 yj , 

where 11 = 506, i? = 306, Li=A0B and L -2 = 60D. Here L\ and L 2 are the periods of 
the perturbations in the x and y directions, respectively, and the wavenumbers of the 
perturbations in DDPFs </> and i/j are (10,2) and (2,5), respectively. 
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Fig. 8.3. Example 3: Continuum glide force and comparison with that of the discrete dislocation 
model for a general dislocation distribution given by Eq. (8.3). The force unit is fib. 


Fig. 8.3 shows the values of the continuum glide force calculated by Eq. (6.12) 
and comparisons with the results obtained using the discrete dislocation model. It 
can be seen that the glide force profile calculated by the continuum model including 
both the long-range and short-range interactions (in Fig. 8.3(a)) excellently keeps the 
overall features of the glide force profile calculated by the discrete dislocation model 
(in Fig. 8.3(c)), whereas the continuum long-range glide force alone (in Fig. 8.3(b)) 
loses too much detailed information compared with the discrete model (in Fig. 8.3(c)). 
Moreover, as shown in Fig. 8.3(d) and (e), the full continuum force successfully reduces 
the maximum error of the continuum long-range force by half, although the continuum 
short-range terms are derived only from special distributions of dislocations. 
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8.2. Dynamics simulations 

In this subsection, we present some simulation results for the dynamics of the dislo- 
cation structures and compare the results with those of the discrete model. We consider 
the dislocation distribution of Case 1 in Sec. 5. In this case, the continuum model is 
given by Eq. (6.13). We fix the uniform active slip plane spacing if> y = D = 50b. 

The initial state of the evolution is a dislocation wall system of N = 40 dislocation 
walls with average spacing B = 30 b. The left half of these dislocation walls consist of 
dislocations with direction in the +z direction, and the right half consist of dislocations 
with direction in the — z direction. Initially, these dislocation walls have equal spacing. 
An initial profile of the DDPF cj> is shown by the blue curve in Fig. 8.4. We use periodic 
boundary condition in the simulations. The dislocation walls at the two ends of the 
simulation domain are fixed. We evolve the dislocation system under applied shear 
stress cr° y = —0.0009 fib and —0.09 fib. 



x/B 


Fig. 8.4. Evolution of the dislocation walls system (represented by the evolution of the DDPF 
ft) and equilibrium locations of dislocation walls (dots on the x-axis) under applied shear stress tj® y = 
—0.0009. The blue curve is the initial profile of (j>(x), and the red curve is the profile of cft(x) of the 
final, equilibrium state. 


Evolution of the dislocation walls system represented by the DDPF <j> under applied 
shear stress <r° y = —0.0009 is shown in Fig. 8.4. It can be seen that during the evolution, 
some opposite-direction dislocation wall pairs initially in the middle annihilate, and the 
remaining dislocation walls are piled-up at the two ends of the domain. Finally, an 
equilibrium state is reached, in which the +z dislocation walls are piled up at the left 
end of the domain and the — z dislocation walls are piled up at the right end of the 
domain. 

The obtained equilibrium dislocation wall distributions under applied shear stress 
a xy = —0.0009 fib and —0.09^6 and comparisons with the results obtained by discrete 
dislocation model are shown in Fig. 8.5. In both cases, the simulation results using 
the continuum model agree excellently with the results of discrete dislocation model for 
these pile-ups of dislocation walls, even though the dislocation wall densities are high 
in the pile-ups and vanishes in the middle of the domain. 

9. Conclusions 

In this study, we have considered systems of parallel straight dislocation walls and 
have identified four cases of these dislocation structures where the continuum long-range 
force vanish but the Peach-Koehler force from the discrete dislocation model does not. 
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(a) Wall locations for a® y = — 0.0009/rfe (b) Wall density for cr® y = — 0.0009/rfe 



Fig. 8.5. Equilibrium dislocation wall pile-ups calculated using our continuum model and compar- 
isons with the results of the discrete dislocation model under the applied shear stress (7? ry = — 0.0009/xb 
(images (a) and (h ) ) and cr? ry = —0.009/d) (images (c) and (d)). The profiles of the DDPF 0 and 
locations of the dislocation walls in the equilibrium states are shown in (a) for (j® y = — 0.0009/rb and 

(c) for cr® y = - 0.009/d). The black dots on the x-axis indicate the locations of the dislocation walls, 
and the blue dots show the corresponding values of <f> in the continuum model. The densities of the 
dislocation walls (given by <f>'(x)) in the equilibrium states are shown in (b) for = —0.0009 pb and 

(d) for a° y = —0.009 pb. 


We have developed continuum descriptions for the short-range dislocation interactions 
for these four cases by using asymptotic analysis. The obtained continuum short-range 
interaction formulas are incorporated in the continuum model for dislocation dynam- 
ics based on a pair of dislocation density potential functions that represent continuous 
distributions of dislocations. This derived continuum model is able to describe the 
anisotropic dislocation interaction and motion. It has been shown that after incor- 
porating these short-range interaction terms, the continuum model is able to provide 
strong stabilizing effect as does by the discrete dislocation dynamics model. Since these 
short-range interaction terms are in the form of second order partial derivatives of the 
DDPFs <p and ip, they also serve as regularization terms in the evolution equations of 
<p and ip. The derived continuum model is validated by comparisons with the discrete 
dislocation dynamical simulation results. 
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